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Abstract: This paper presents a complete procedure for sensor compatibility correction of 
a fixed- wing Unmanned Air Vehicle (UAV). The sensors consist of a differential air 
pressure transducer for airspeed measurement, two airdata vanes installed on an airdata 
probe for angle of attack (AoA) and angle of sideslip (AoS) measurement, and an Attitude 
and Heading Reference System (AHRS) that provides attitude angles, angular rates, and 
acceleration. The procedure is mainly based on a two pass algorithm called the 
Rauch-Tung-Striebel (RTS) smoother, which consists of a forward pass Extended Kalman 
Filter (EKF) and a backward recursion smoother. On top of that, this paper proposes the 
implementation of the Wiener Type Filter prior to the RTS in order to avoid the 
complicated process noise covariance matrix estimation. Furthermore, an easy to 
implement airdata measurement noise variance estimation method is introduced. The 
method estimates the airdata and subsequently the noise variances using the ground speed 
and ascent rate provided by the Global Positioning System (GPS). It incorporates the idea 
of data regionality by assuming that some sort of statistical relation exists between nearby 
data points. Root mean square deviation (RMSD) is being employed to justify the sensor 
compatibility. The result shows that the presented procedure is easy to implement and it 
improves the UAV sensor data compatibility significantly. 
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1. Introduction 

Air vehicle sensor compatibility is a specific field of study dedicated to the estimation of onboard 
sensor errors. The sensor errors are usually described as the scale factors and the biases. The field is 
often discussed in the context of flight path reconstruction. The term "compatibility" refers to the 
consistency of the measured data with the force equations that govern the airdata (airspeed, AoA, and 
AoS), and the kinematic equations that govern the attitude angles (roll, pitch, yaw). If flight path 
reconstruction is to be discussed jointly, the navigational equations will be included. Sensor 
compatibility is often the first step of the so-called "two-step method" of aircraft system identification. 
The first step is the sensor data correction/estimation; the second part is the parameter estimation. 

Research on the topic of aircraft sensor compatibility started as early as the 1970s. In 1976, 
Jonker [1] pioneered the research in the area by investigating the applicability of the Kalman Filter for 
flight path reconstruction and estimation of instrumentation biases. A year later, the National Aeronautics 
and Space Administration (NASA) published a report on the measured aircraft responses compatibility 
check. As discussed in the report, Klien and Schiess [2] implemented the EKF with a fixed-point 
smoother to accomplish the task. In 1984, similar work was again published by NASA. This time, 
Whitmore et al. [3] calibrated the airdata error using the Linearized Kalman Filter (LKF). Up to that 
point, The Kalman Filter and EKF were proven to be sufficient in estimating the air vehicle sensor 
errors. However, the implementations of Kalman Filter and EKF are not easy without a priori 
knowledge of the measurement noise and process noise and this noise information is often not easily 
available. In response to that issue, Chu et al. [4] proposed the implementation of the modified 
recursive Maximum Likelihood adaptive filter in 1995, claiming that the adaptive filter is more robust 
than the EKF. Mendonga et al. [5] also proposed an adaptive method for estimating the noise 
properties in 2007. All these efforts were carried out using manned aircraft. 

In the age of the UAV, due to the advances in discrete time system identification techniques and 
control theories, identifying an accurate model that is coherent with the force equations and kinematic 
equations seems to be unnecessary. However, discrete time system identification techniques are unable 
to provide a good insight of an aircraft's dynamics because the discrete time models are usually 
physically meaningless, no matter how well they agree with the real flight data. In addition to that, 
UAVs are usually smaller than manned aircrafts, and fly at lower velocities, and thus are very sensitive 
to the atmospheric turbulence which will result in a higher noise ratio in the flight data. Due to that, 
sensor compatibility correction is more important than ever if the dynamic characteristics of an UAV 
are to be investigated. 

Motivated by the discussion above, this paper attempts to derive a procedure in order to perform 
compatibility correction on the flight data of the SP-80 UAV. The SP-80 UAV [6-9] is an UAV 
research testbed developed and operated under the Spoonbill UAV Project by the Remotely Piloted 
Vehicle and Micro-Satellite Laboratory (RMRL) of the National Cheng Kung University (NCKU), 
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Taiwan. The Spoonbill UAV Project was established to fulfill educational research purposes of which 
it has brought forth the successful implementation of subspace and prediction-error method 
identification algorithms in the UAV system identification [7], implementation of the Linear Quadratic 
Gaussian (LQG) controller to the autopilot system [8] and a successful cross-sea automatic flight 
demonstration that covered 90.2 km [6]. As part of the continuing researches on the SP-80 UAV, the 
work presented in this paper intends to bridge the gap between the controller implementation and the 
understanding of the dynamics of the UAV. 

The following sections in this paper are organized as follows: first, Section 2 introduces the Spoonbill 
UAV system. In particular, the sensor system is described in detail. Section 3 reviews the force equations 
and kinematic equations. This is followed by Sections 4 and 5 which briefly elaborate the Wiener Type 
Filter and the RTS smoother. Section 6 presents the flight test data and initial data reconstruction and the 
next section explains the sensor compatibility correction procedure in detail. Then, the result of 
correction is discussed in Section 8. Finally, Section 9 draws the conclusions for the paper. 

2. The Hardware Architecture 

2.7. The Spoonbill UAV System 

The Spoonbill UAV system used in this work is designed and operated by the RMRL for 
educational research purposes. The air vehicle used is designated as the SP-80, which has a 
high-wing, twin-boom, pusher engine configuration and the power plant is an 80 c.c. two stroke 
gasoline engine. Figure 1 depicts the air vehicle. The 3.50 m wingspan air vehicle is capable of 
carrying 3 kg of fuel and it has a maximum endurance of an hour once fully refueled. The fuselage is 
designed such that it can carry various sensors, an onboard computer, telemetry devices, and a 
gimbaled camera system. So far, the air vehicle is capable of performing automatic flight using a flight 
controller based on the LQG algorithm, and it automatically flew the longest distance of 90.2 km on 20 
October 2009 [6]. 



2.2. The Onboard System 

Figure 2 shows the onboard avionics system and Figure 3 illustrates the Spoonbill UAV system 
architecture. The onboard computer acts as the central processing unit. It is a PC/104 embedded 
computer produced by Advantech Co., Ltd. The onboard computer is operating under a 
Windows-XP embedded operating system and the flight control computer program is written using 



Figure 1. The SP-80 UAV. 
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Borland C++ Builder. The onboard computer collects data from the peripheral sensors through RS-232 
protocol via several communication ports (COM Port). In addition to that, the onboard computer also 
performs bi-directional communication with the Servo Management Board (SMB), and the wireless 
transceiver. The final COM Port is occupied by the rotation platform of the gimbaled camera system. 
The flight control program operates under various modes depends on the system operator and the flight 
test objectives. The basic functions of the flight control program are: acquires data from the peripheral 
sensors, executes the automatic flight control algorithm, and records selected flight data. All of these 
operations are executed at the rate of 20 Hz. 

Figure 2. The onboard system. 




Figure 3. Spoonbill UAV System architecture. 
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Before any discussion on the peripheral sensors can commence, it is worth mentioning that the GPS 
receiver, the AHRS, the wireless transceiver, and the video camera system are all off-the-shelf 
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products, while the Sensor Integration Board (SIB), and the SMB were designed and fabricated 
in-house by the RMRL. The SMB is a crucial component of an autopilot system for it decodes the 
control signal from the pilot in Pulse Width Modulation (PWM) form during the manual flight mode, 
and it encodes PWM signal to control the air vehicle control surface actuators (servo motors) 
according to the command given by the onboard computer. The SMB is capable of decoding and 
encoding nine channels of PWM signals at a 20 Hz update rate. The pilot can easily cede flight control 
authority to the onboard computer by flipping a switch on the remote control radio. Regaining control 
of the air vehicle from the onboard computer is just as easy. 

The peripheral sensors are divided into three major groups. The first group is the GPS. The GPS 
receiver is the Novatel OEMV-3. Unlike any other GPS receivers built for automobile navigation that 
normally operating at a maximum of 5 Hz, The Novatel OEM-V3 is capable of providing global 
position and ground speed at the rate of 20 Hz and is thus consistent with the sampling rate of other 
onboard sensors. The second group is the AHRS. It is a Crossbow AHRS440. The AHRS measures 
3-axis accelerations, angular rates, and attitude angles. The attitude angle measurement is augmented 
by GPS information from the Novatel OEM-V3. The final group is the SIB. The SIB integrates all 
other sensors using microcontroller units (MCU). It is designed and fabricated in-house by the RMRL. 
The MCUs on the SIB are responsible to collect data from the corresponding sensors, perform signal 
conditioning on the data (such as data filtering), and pass the data to the onboard computer via RS232 
protocol. Sensors included in the SIB are the differential air pressure sensor for airspeed measurement, 
absolute air pressure sensor for barometric altitude measurement, engine speed tachometer, AoA and 
AoS, and the touchdown sensors. Figure 4 illustrates the SIB architecture. 

Figure 4. Sensor Integration Board architecture. 
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The SIB consists of two MCUs: the Master MCU and the Slave MCU. The Master MCU is the core 
of the SIB. It controls the data flow within the SIB and sends data to the onboard computer. There are 
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two sensors built on the SIB: the SenSpecial SCPB-MB0/10D100i2c32667R5 digital differential 
pressure sensor, and the SenSpecial SCPB-mmHg525/825A100i2c32667R5 digital absolute pressure 
sensor. Both pressure sensors are connected to the air pressure tubes from the SpaceAge Control 
101100-02 airdata probe. The airdata probe is mounted on the starboard wing, as shown in Figure 2. 
The Master MCU communicates with the pressure sensors via I C protocol, and calculates the airspeed 
and barometric altitude using the pressure readings and the following equations. 



\ P b J 



g 0 UM 



-1 



V = 



P 



p 

RT 



T = T b +L b h 

Where g 0 = 9.80665m/ s 2 is the gravitational acceleration at sea level 
h is the geopotential altitude in meter 
h b = 0 is the sea level altitude 

L b = -6.5 x 10~ 3 K/m is the temperature lapse rate from sea level to altitude 1 1 ,000 m 

M = 0.0289644 kg/mol is the molar mass of Earth's air 

P is the static pressure in Newton per square meter (Pascal) 

P 0 is the total pressure in Newton per square meter (Pascal) 

P b = 1.01325 x 10 5 n/ m 2 is the sea level static pressure 

R = 287 J I (kg - K) is the specific gas constant 

/?* = 8.314327V • ml (mol - is the universal gas constant 

T is the temperature 

T b = 288A6K is the sea level temperature 

V is the airspeed in meter per second 

p is the air density in kilogram per cubic meter 



(1) 

(2) 

(3) 
(4) 



Equation 1 is called the barometric formula introduced by NASA in 1976 [10]. It has since become 
the standard equation for barometric altitude calculation. Equation 2 is derived from the Bernoulli's 
equation [11]. From Equation 2, the airspeed of the air vehicle can be determined once the total 
pressure, static pressure and air density are known. Due to the absence of an air density sensor, the air 
density is estimated using Equation 3. The temperature needed in Equation 3 is calculated by one of 
the standard atmosphere formulae shown here as Equation 4 [10]. Equation 4 is only valid from sea 
level to altitude 11,000 m. The absolute pressure sensor measures P , while the differential pressure 
sensor measures (P 0 - P) . 

The airdata probe is also built with two vanes, one for the AoA measurement, and another for the 
AoS. The vanes are basically two delicate potentiometers. The SIB provides the operating voltage to 
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the vanes and senses the vanes' movement through the Master MCU 12 bit analog-to-digital 
converters (A/D). 

The onboard avionics are powered by a high capacity Lithium-Polymer (Li-Po) battery and the air 
vehicle control surface actuators (servo motors) are powered by high capacity Nickel-Cadmium 
(Ni-Cd) batteries. The Slave MCU reads the battery voltages via a built-in 10 bit A/D. In addition, it 
reads the pulse signal from the touchdown switches mounted on both main wheels of the air vehicle; 
such to determine the instant the air vehicle lifts off the ground and touches down. The most important 
function of the Slave MCU is to calculate the engine rotation speed. The opto interrupter mounted near 
the engine shaft generates a pulse per engine revolution. The Capture/Compare/PWM (CCP) interrupt 
of the Slave MCU detects the pulse and calculates the time taken between two pulses before converts it 
to revolution per minute (rpm). The Slave MCU sends the data string to the Master MCU upon Master 
MCU'sI 2 C request. 

3. Equations of Motions 

The air vehicle translation and rotation dynamics are governed by two sets of ordinary differential 
equations (ODE), namely the force equations and the kinematic equations. Each set of these equations 
consists of three equations, and each equation governs the motion of one axis. The detailed derivation 
of the equations has been described well by Nelson [12]. 

3.1. Force Equations 

The force equations as shown in Equation 5 describe the variation of the three body axis velocities. 
The ODEs are functions of the velocities, accelerations, angular rates, and gravity components: 



where: a x is the acceleration along the body X-axis 
a y is the acceleration along the body Y-axis 
a z is the acceleration along the body Z-axis 

_9 

g = 9.8 lms is the gravitational acceleration 
p is the roll rate 
q is the pitch rate 
r is the yaw rate 

u is the velocity along the body X-axis 
v is the velocity along the body Y-axis 
w is the velocity along the body Z-axis 



u = rv - qw + a x - g sin 0 



(5a) 



v = pw -ru + a y + g cos 6 sin $ 
w = qu- pv + a z + g cos dcos <fi 



(5b) 



(5c) 
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It is also possible to describe the force equation in terms of the variation of airspeed, AoA, and AoS 
The airspeed, AoA, and AoS are related to the body axis velocities as in Equation 6. Thus, the force 
equations can be transformed into Equation 7 that describe the variation of airspeed, AoA, and AoS: 

V = ^lu 2 +v 2 +w 2 



a = tan 



P = sin" 1 



V = (a x cos a + a z sin a) cos J3 + a sin ft + 

g (cos 9 cos (j) sin a cos /? + cos#sin^sin/?-sin<9cosacos /?) 

' \a z cosa-a x sina + g (cos # cos ^ cos a + sin6 ) sina)} 



a = 



V cos p 
q-tanp(p cos a + r sin a ) 

a cosP -(a x cosa + a z sin a) sin P + 

g[cos^sin^cos y^ + (sin^cosa - cos ^ cos ^ sin a) sin /?] 



+ 



+ psina-rcosa 



(6a) 
(6b) 

(6c) 
(7a) 

(7b) 
(7c) 



where: V is the airspeed 
a is the AoA 
P is the AoS 

As a matter of fact, even though Equation 7 is much more complicated compared to Equation 5, it is 
more intuitive to use Equation 7 rather than Equation 5 because the airspeed, AoA, and AoS are 
normally measured by the onboard instruments instead of using the body axis velocities. 

3.2. Kinematic Equations 

The kinematic equations as shown in Equation 8 describe the variation of the three attitude angles, 
namely the roll angle, pitch angle, and yaw angle. The ODEs are functions of the attitude angles, and 
angular rates: 

/? + tan#(gsin^ + rcos^) 



9 = qcos<p-rsin<p 
q sin (j) + r cos <f> 



COS0 



(8a) 
(8b) 

(8c) 



where: $ is the roll angle 

9 is the pitch angle 
(j) is the yaw angle 
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4. Wiener Type Filter 

The Wiener Type Filter [13] is an optimal filter in the frequency domain. Equation 9(a) indicates 
that the measurement is a summation of the filtered/true signal and noise. In order to simplify the 
transformation to the frequency domain, the measurement is modified using Equation 9(b). 
Equation 9(b) is an odd function with its end point discontinuity removed. Such function allows the 
Fourier sine series expansion. The Fourier sine series coefficients are given by Equation 9(c): 

*k =y k + n k> k = l 9 2,...,N (9a) 



N-l 



b,=- 



N- 



A 1=2 



sin 



lx(k-l) 
N-l 



, £ = 1,2,. ..,7V 



, / = 1,2,...,7V-1 



(9b) 



(9c) 



where: b t is the Fourier sine series coefficient at I th frequency index 



g k is the modified measurement at k th time step 
N is the data size 



n k is the noise at k th time step 



y k is the filtered/true at k th time step 
z k is the measurement at k th time step 

The next step is to define the true signal model and noise model. It is important to know that the 
filter consists of weighting indices range from 0 to 1, where 0 indicates total rejection and 1 indicates 
otherwise. Thus, it is ideal when the filter equals to 0.5 as the frequency index then corresponds to the 
cutoff frequency. Plus, as discussed by Morelli, it is best to describe the true signal model being 
proportional to /" 3 . Given the form of the Wiener Type Filter as in Equation 10(a), yields the true 
signal model and noise model as shown in Equation 10(b,c). Equation 10(d) shows how the frequency 
index being related to the frequency. Similarly, Equation 10(e) determines the frequency index 
corresponds to the desired cutoff frequency. The filter is applied in Equation 10(f) to obtain the filtered 
modified measurement. Finally, Equation 10(g) reconstructs the filtered measurement: 



O, =■ 



, 0<1<N-1 



y, = Cr 



N t =l 



//=■ 



>s,k 



2(N-l)At 
/ c =2/ l (iV-l)Af,if f t <f c4 <f M 
^ 'lx(k-l) 



sin 



i=i 



N-l 



, k = l2,...,N 



(10a) 

(10b) 
(10c) 

(lOd) 

(lOe) 

(lOf) 
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y k =8s, k +z 1 +(k-l) 



N-l 



(10g 

) 



where: f x is the frequency at / frequency index 
f c d is the desired cutoff frequency 
g s k is the filtered modified measurement at k th time step 
/ c is the frequency index corresponds to the cutoff frequency 
N l is the noise model in discrete frequency domain at I th frequency index 
Y t is the true signal model in discrete frequency domain at I th frequency index 
O t is the filter at I th frequency index 

5. Rauch-Tung-Striebel Smoother 

The RTS Smoother [14-16] is a two pass algorithm. The forward pass is a standard EKF [15] while 
the backward recursion is introduced to reduce the inherent bias in the EKF estimates. 

5.1. Forward Pass (Extended Kalman Filter) 

The EKF is literally one of the most widely used state estimation methods due to its simplicity. The 
EKF is derived from the standard Kalman Filter. The Kalman Filter is derived based on linear systems, 
but the EKF works for nonlinear systems by performing linearization about the current estimates. 

Over the years, different forms of EKF have been proposed, namely the continuous time EKF, 
discrete time EKF, and the continuous-discrete EKF. Since the force equations and kinematic 
equations are continuous time equations and the flight test measurements are discrete, the 
continuous-discrete EKF is used for the work presented in this paper. The continuous-discrete EKF is 
presented here as Equations 11-12. Often, the difficulty of implementing the EKF is mainly due to the 
lack of a priori knowledge of the noise covariance. 

Given a nonlinear continuous time dynamic system with discrete time measurement, of which the 
measurement and process noise are Gaussian distributed: 



h(x k ) is the measurement equations 

k is the number of time step 

Q is the process noise covariance matrix 

R is the measurement noise covariance matrix 

u is the input vector 

v is the Gaussian distributed measurement noise vector 
w is the Gaussian distributed process noise vector 
x is the state vector 
z is the measurement vector 




(lib) 



(lla) 
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The prediction phase: 



t+At p^r 

i(f)= \ f(*(t),n(t))dt,V(t) = CJ 



The updating phase: 



where: P(f t ) = P, 
P('*-i) = P| 



k\k-\ 

k-\\k-\ 
v k\k-l 



^k-\\k-\ 



dx 



x(t),u(t) 



p(0=F(r)p(0+p(0 F W r + Q(0 



K * - P *i*-i H / ( H * P *l*-l H / + R * ) 



_ z k z k 



X k\k X k\k-l ~^~^kfk 



P £l£ ~ ^*k\k-l ^k^k^klk- 



F is the Jacobian of the nonlinear ordinary differential equations 

H is the Jacobian of the measurement equations 

K is the gain matrix 

P is the state error covariance matrix 

x is the estimated state vector 

y is the measurement residual 

z is the estimated measurement 



(12a) 
(12b) 
(12c) 

(12d) 
(12e) 
(12f) 
(12g) 



Equation 12(b) is a continuous time Riccati equation. In order to solve the equation, the discrete 
state transition matrix can be determined using Equation 12(h). With the state transition matrix 
available, Equation 12(i) provides the solution of the state error covariance matrix. Note that 
Equation 12(h) is applied under the assumption that F is a constant from t k _ x to t k : 

(D H ^* (12h) 

= ®k-l?k-M-l®ll + Q*-l (12i) 

where <l> is the state transition matrix 

It is important to know that the EKF needs predefined P 0I0 , R , Q , and x 0l0 for initialization. It is 
acceptable to assign P 0I0 arbitrarily, however, the values must be large enough to allow good tracking 
of the parameters. If the states are measured, x 0l0 can be specified by taking the average of the first 
few data points. 
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5.2. Backward Recursion 



The backward recursion works backwards in time. The combination of a forward pass estimator 
(EKF) and a backward recursion is considered to have utilized all available information [15]. Thus the 
scheme is capable of performing better estimations than the EKF. Equation 13 summarizes the 
backward recursion of the RTS smoother: 

A,=P^/P, + m _1 (13a) 
Kn =Ku+ K(K + i\ n -K + i\k)> k = N-l... 9 0 (13b) 

^*k\n - ^*k\k + ^k [j*k+l\n ~ ^k+l\k ^^k (13C) 

where: A is the smoother gain matrix 
N is the final time step 

P kln is the corresponding state error covariance matrix 
x kln is the smoothed states of k th time step 



6. Flight Tests 

6.1. The Flight Maneuver 

The flight test data were collected through a series of flight tests. In the flight data acquisition 
process, it is desirable to maneuver the air vehicle such that the flight data contain sufficient 
information to represent the dynamics of the air vehicle. The maneuvers chosen to fulfill the objective 
of this work were a multiple-input design. Due to its simplicity, time-skewed doublet inputs like those 
used in 2008 by Lee et al. [7,8] are always a favorable choice, but in order to effectively excite the 
sufficient dynamic information, including the coupling of the longitudinal and lateral motion, the 
orthogonal square- wave inputs [17] were applied. The orthogonal square- wave inputs are depicted in 
Figure 5. 

As shown in Figure 5, at any instant, two control surfaces are deflected simultaneously. One is 
deflected at half the rate of another. The simultaneous deflections are meant to excite couple dynamics, 
while the separate deflection rates are meant to distinguish the corresponding control surfaces' effects 
on the air vehicle motion. 

The initial condition of the designed flight maneuver has to be trimmed flight condition. Due to the 
complexity of the maneuver and the difficulty for the ground pilot to ensure a trimmed flight initial 
condition, the maneuver was executed by the onboard computer. The ground pilot flew the SP-80 
UAV to a desired altitude and remained straight and level flight before switching the control authority 
to the onboard computer. Once switched, the onboard computer maintained the air vehicle in trimmed 
straight and level flight for 5 seconds. After that, the designed maneuver was executed and the ground 
pilot switched it back to manual flight no less than 5 seconds after the maneuver. The desired altitude 
was about 300 m above ground. Since the air vehicle remained at open-loop control during the 
maneuver, the 300 m altitude gave the ground pilot enough time to react to any unexpected outcome 
while the air vehicle remained in the pilot's eyesight. 
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Figure 5. Orthogonal square- wave inputs. 
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It is desirable to perform the flight test under absolutely windless conditions, but this is utterly 
impossible. Thus, the flight tests were conducted in nearly windless environment and the maneuvers 
was executed in both head wind and tail wind directions to minimize the crosswind effect. 

6.2. The Flight Data and Reconstruction 

The flight data reconstruction is accomplished by implementing the 4th order Runge-Kutta method. 
The accelerations and angular rates are the inputs; the airdata, and attitude angles are the outputs. The 
4th order Runge-Kutta method is summarized as Equation 14: 

k i =h-f(x(t) 9 u(t)) 



k 2 =h-f 



k 3 =h-f 



*(<)4 



x(t) + ^ 




(14) 



k 4 =h-f([x(t) + k 3 ],u(t)) 

x(t + h) = x(t) + k l +2k 2+ 2k l+K 



where h is the time step size 

Multiple sets of flight data were collected and each data set is designated in this paper in the 
following format: ftFFsrtSSrunRR, where FF is the flight test number of the Spoonbill UAV project, 
SS is the sortie number, and RR is the number of maneuver in that particular sortie. This paper 
presents two data sets as examples: ft74srt02run01 and ft74srt02run05. Data of ft74srt02run01 is 
depicted in Figure 6, and Figure 7 corresponds to ft74srt02run05. Both data sets show similar trends. 
The reconstructed attitude angles demonstrate very good compatibility while the airdata compatibility 
is poor. 

From Equation 7, it is legitimate to infer that the airspeed, AoA, and AoS variations are dominated 
by X-axis acceleration, Z-axis acceleration, and Y-axis acceleration, respectively. The reconstructed 
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airspeed is over estimated in both data sets. It is obvious that the reconstructed airspeeds are increasing 
at a higher rate compared to the measured airspeed, which suggest the presence of a bias on the X-axis 
acceleration. Fluctuations are observed on the measured airspeeds at the 9th to 12th second, but not on 
the reconstructed airspeeds. That could be the result of instrumentation errors of the airdata probe due 
to rapid yawing motion of the air vehicle. The AoA and AoS demonstrate slightly better 
compatibilities compared to the airspeed. The reconstructed data seem to follow the trend of the 
measured data. This suggests that the biases of Y-axis acceleration and Z-axis acceleration are less 
significant. Thus, the instrumentation error, measurement noises, and atmospheric turbulent are the 
major causes of the AoA and AoS incompatibilities. 

Similar inferences were made on Equation 8. The roll angle, pitch angle, and yaw angle variations are 
dominated by the roll rate, pitch rate, and yaw rate, respectively. It is well understood that the attitude 
angles are estimated through complex algorithms using the angular rates, and accelerations [18]. The 
principle idea is that the attitude angles are estimated through integration of the angular rates, and the 
accelerations prevent the drifting of the estimation. That explains the high compatibilities of the 
attitude angles. Since the attitude angles estimation is augmented by GPS data [19,20], it is more likely 
that the attitude angles are relatively accurate and noises of the angular rates are the reason of the 
slight incompatibilities. 

Figure 6. (a) Accelerations and angular rates of ft? 7 4 srt02runOL (b) Airdata and attitude 
angles of ft74srt02runOL 
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Figure 6. Cont. 
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Figure 7. (a) Accelerations and angular rates of ft74srt02run05. (b) Airdata and attitude 
angles of ft74srt02run05. 
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Figure 7. Cont. 
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7. Sensor Data Compatibility Correction 

7.1. The Compatibility Correction Procedure 

The sensor compatibility correction procedure as illustrated in Figure 8 was designed under several 
assumptions: 

• The AHRS body axes and the air vehicle body axes coincide. 

• The airdata probe is located near to the center of gravity such that position correction is 
unnecessary. 

• Since the AHRS is not exposed to the environment, the AHRS readings (accelerations, 
angular rates, and attitude angles) are less affected by atmospheric turbulence, and the 
noises (not the biases) are high frequency noises that can be effectively filtered by a 
frequency filter. 

• Due to the high compatibility of the attitude angles, the biases of the angular rates are 
negligible. 

• The noises or errors of different variables are not inter-correlated. 

• The noises distributions are Gaussian. 

As seen in the right portion of Figure 8, the flight data undergo reconstruction using the 4th order 
Runge-Kutta Method. The reconstructed results of ft74srt02run01 and ft74srt02run05 are presented in 
the previous section as Figure 6 and Figure 7. The RMSD shown here as Equation 15 is then 
implemented to quantify the quality of compatibility. The smaller the RMSD is, the better the 
compatibility is. 
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Figure 8. The compatibility correction procedure. 
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The RMSD between the measured and reconstructed airdata and attitude angles are calculated. This 
is the RMSD before the correction: 



RMSD(y 19 y 2 ) = 



j=i U^j 



n d 



where: n d is the data size 



Yj is the first data 

y 2 is the second data 

Y u is the i th element of the first data 

y 2 i is the i th element of the second data 



The left portion of Figure 8 illustrates the correction process. The RTS smoother is the core of the 
whole process. The flight data (airspeed, AoA, AoS, attitude angles, angular rates, and accelerations) 
is first passed through a Wiener Type Filter to suppress the high frequency noises of the flight data. It 
then followed by the airdata measurement noise estimation that will be discussed in depth in 
Section 7.3. RTS is implemented after that to correct the flight data such that the data is compatible to 
the force equations and kinematic equations. The corrected accelerations and angular rates are then 
used as inputs to perform airdata and attitude angle reconstruction. If successful, the corrected and the 
reconstructed airdata and attitude angles should show a high degree of similarity. After that, the 
RSMD between the corrected and the reconstructed outputs (airdata and attitude angles) are calculated. 
This is the RMSD after the correction. Finally the RMSDs before and after the correction are 
compared to demonstrate the sensor data compatibility improvement. For the work presented in this 
paper, the RTS was executed with 0.05 s sampling time, which is consistent with the 20 Hz sampling 
rate of the onboard system. 
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7.2. RTS Problem Formulation 



Under the assumptions discussed in the previous section, the RTS problem formulation is as shown 
in Equations 16-18. Equation 16 shows the relation between the measured accelerations and the 
corresponding biases. As in Equation 17, the biases are formulated as the state variables such that the 
RTS smoother could estimate them. The airdata and the attitude angles are part of the state variables 
and they are also the measurement. The accelerations and angular rates are the inputs. While the force 
equations and kinematic equations serve as the ordinary differential equations of corresponding states, 
the acceleration biases are assumed constant as shown in Equation 18: 

= <2 + b 



(16a) 



a y M = a y+K 



(16b) 



= a+b 



where: a v is the measured X-axis acceleration 



X M 



a is the measured Y-axis acceleration 



a v is the measured Z-axis acceleration 

b a is the X-axis acceleration measurement bias 

b is the Y-axis acceleration measurement bias 

b„ is the Z-axis acceleration measurement bias 



(16c) 



x = [V a J3 (j> 0 y/ bf,b = 



z = [V a P (j) 0 y/f 



u = [a x 



yM 



pgr 



where b is the acceleration measurement bias vector: 

b = (h 



'3x1 



(17a) 
(17b) 
(17c) 

(18) 



73. Measurement Noise & Process Noise Covariance Matrix Estimation 

It was assumed that the noises are not inter-correlated in order to simplify the noise covariance 
matrices estimation. Hence, the task is to estimate the diagonal terms of the matrices, namely the 
variances. The measurement noise variances of the airdata are estimated from the speeds measured by 
the GPS. Unlike the airdata, the GPS speed measurements are less noisy. That is because the GPS 
doesn't measure speed by interacting with the relative airflow like the airdata probe does. In addition, 
the atmospheric turbulence is the major contributor to the airdata noise or error. One might argue that 
it is impossible to accurately estimate the airdata measurement noise variances from the GPS speed 
measurements because the airdata are relative to the wind, but the GPS measurements are relative to 
the ground. But it is later demonstrated that the GPS speed measurements are sufficient to fulfill the 
purpose. 
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The GPS speed measurement would have to go through a series of coordinate transformations 
before any noise variances estimation is possible. Equation 19 transforms the horizontal ground speed 
and ascent rate measured by the GPS to the North-East-Down (NED) local horizontal coordinate frame. 
The NED speed components are then transformed to the air vehicle body-axis using Equation 20. After 
that, the body-axis speeds are substituted into Equation 21 to obtain estimations of total speed, AoA, 
and AoS. It is essential to take note that these estimations are not meant to be accurate, but since they 
are far much less sensitive to atmospheric turbulent, these estimations provide a good baseline 
reference for the airdata measurement noise variances estimation: 



V = V cos w 

N,gps hor,gps t 



y = y sin ur 

E,gps hor,gps t 



V = -V 

D,gps U,gps 



where: V hor gps is the horizontal speed given by GPS 
V E is the ground speed to the East 
V N is the ground speed to the North 
V U gps is the ascent rate given by GPS 
V n m , is the descent rate 

u tgPS 



(19a) 
(19b) 
(19c) 



u 




gps 




V 




gps 




w 









1 0 0 

0 cos <f> sin <f> 
0 -sin^ cos^ 



cos0 0 



0 
sin# 



1 
0 



-sin# 
0 

cos 6 



cosy/ 
-siny 
0 



sin^ 
cos^ 
0 



V, 



N,gps 



V, 



E,gps 



V, 



D,gps 



(20) 



where: u is the body X-axis speed estimated from GPS measured speeds 
v is the body Y-axis speed estimated from GPS measured speeds 
w is the body Z-axis speed estimated from GPS measured speeds 



gps 



2 2 2 

+ v +w 

gps gps gps 



V U sps J 



gps 

, V , 



(21a) 
(21b) 

(21c) 



where: V is the GPS total speed 

gps -L 

a onv is the GPS estimated AoA 

gps 

B m is the GPS estimated AoS 

/ gps 

The airdata estimations given by Equation 21 serve as reference for the airdata measurement noise 
estimations. The most intuitive way of estimating the noise is determining the residual between the 
airdata estimations with the airdata measurements. However, the airdata estimations come with 
inherent biases due to the difference between airspeed and ground speed. Instead, the measurement 
noise estimations are derived based on the deviation of the airdata at any particular time with their 
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local expected values. The idea is as shown in Equation 22. Take Equation 22(a) as an example, the 
first portion of the equation tells how far is the k th airspeed measurement deviate from its local 



expected value. The local expected value is the mean of airspeed measurements from (k - m) th to (k + 



m) th time steps. Similarly, in the second portion of the equation, the deviation of the k m GPS total 
speed from its expected value is calculated. The difference between these two deviations gives the 
estimated noise of the k th airspeed measurement. One can control the local expected values by 
changing m. If the air vehicle is highly dynamic, a smaller m gives more legitimate local expected 
values: 



,th 



S V,k 



S a,k 



1 2m+ \ 1 2m+ l 

V v. - V V V 

2m + lit 'I I ~* 2m + l£ m t 



J 2m+l J 2m+l 

~2m + l " ^ _ 2m + l,£ 



1 2m+l 

— £ Pi 

2m + l i=k _ m 



0, 



gps,k 



1 2m+l 

— — y ^ 

2m + l /= £^. 



(22a) 



(22b) 



(22c) 



where: e is the estimated noise/error 
m is the design parameter 
Once the airdata measurement noises are estimated using Equation 22, the noise variances can be 
calculated. Using the same idea of data regionality, the noise variances are calculated by taking the 
nearby data points into account. The equations are presented as Equation 23. Taking Equation 23(a) as 
an example, it is a standard variance equation; the k th airspeed measurement noise variance is 
calculated by accounting a sample with sample size of (2 m + 1). The larger the m is, the larger the 
sample size is, and thus the less intense the variation of the noise variance is. Figures 9 and 10 depict 
the estimated airdata measurement noise variance of ft74srt02run01 and ft74srt02run05 at m = 5: 

j 2m+l 2 j 2m+l 



J 2m+\ ^ 

2m + l^_ m v 

j 2m+l 2 

<T ^ = 2m^T, = ?J^"^) ' 



J 2m+l 

= 2mTl£/^') 

J 2m+l 



(23b) 



(23c) 



where: <r is the estimated error variance of airspeed 
<j] is the estimated error variance of AoA 
at is the estimated error variance of AoS 
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Figure 9. Airdata noise variance of ft? 7 4 srt02runOL 
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Figure 10. Airdata noise variance of ft7 4 srt02run0 5. 
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Under the assumptions given in Section 7.1, the measurement noise covariance is a diagonal matrix 
and the attitude angles measurement noise variances are zeros after the implementation of the Wiener 
Type Filter. Hence, the measurement noise covariance matrix is in the following form: 
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0 



0 



0 



0 




0 



3,3 



(24) 



0 



3,3 



0 



3,3 



On the other hand, the process noise is mainly contributed by the accelerations and the angular rates. 
If the second assumption in Section 7.1 is valid, the process noise covariance matrix shall be a zero 
matrix since the data has already been filtered by the Wiener Type Filter. However, a zero process 
noise covariance matrix will lead the RTS smoother to a singularity. Thus, the process noise 
covariance matrix was chosen to be a diagonal matrix with very small diagonal terms as shown in 
Equation 25: 



8. Results and Discussion 

Figures 11 and 12 present the results of the sensor compatibility correction of ft? 7 4 srt02run01 and 
ft74srt02run05, respectively. The figures contain correction results using the EKF (the RTS forward 
pass), the RTS (the RTS forward pass and backward recursion), and the measured data (after Wiener 
Type Filtering at cutoff frequency of 2 Hz). As shown in Figure 11(a) and Figure 12(a), the EKF took 
some time before the acceleration biases converged, and the biases remain constant in the backward 
recursion. This is due to the inherent characteristics of the RTS backward recursion whereby it is 
unable to update the estimations with zero change rates. However, the biases are meant to be constants, 
and as long as the EKF provide converged estimations, the biases are valid and the RTS backward 
recursion helps applying the converged biases to all data points. 

High degree of similarities is demonstrated on the angular rates and attitude angles. Comparing 
Figure 6 and Figure 11, and Figure 7 and Figure 12, the compatibilities of the attitude angles improved 
slightly. The EKF and RTS have done little correction on the attitude angles; the result is mainly due 
to the high frequency noises suppression by the Wiener Type Filter. This is coherent with the 
assumption that the AHRS measurement contains only high frequency noises that can be effectively 
filtered by a frequency filter, and the biases of the angular rates are negligible. Also demonstrated is 
the significant improvement of the airdata compatibility. The reconstructed airspeed tracks the 
corrected airspeed with good consistency. The AoA and AoS also show significant improvement. One 
of the major reasons of the significant improvement of the airspeed compatibility is that the RTS has 
successfully identified a significant bias on the X-axis acceleration as suggested in Section 6.2. The 
improvements on the AoA and AoS compatibilities are due to the noise reduction on the angular rates 
and the biases estimation on the corresponding accelerations. 



Q A =MxlO 



(25) 
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Figure 11. (a) Accelerations and angular rates of ft74srt02run01 before and after 
correction, (b) Airdata and attitude angles of ft? 7 4 srt02run01 before and after correction. 
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Figure 12. (a) Accelerations and angular rates of ft74srt02run05 before and after 
correction, (b) Airdata and attitude angles of ft7 4 srt02run0 5 before and after correction. 
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Tables 1 and 2 demonstrate the quantitative improvement by comparing the RMSDs before and 
after correction. The compatibilities of all the airdata and attitude angles have improved. The airspeed 
has the best improvement with RMSD reduction of as high as 97.76%, which is consistent with the 
significant X-axis acceleration bias estimated in both ft74srt02run01 and ft74srt02run05. Although the 
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improvements of the attitude angles compatibilities are not apparent in Figures 11 and 12, all of the 
attitude angles have shown at least 47.43% of RMSD reduction. 



Table 1. Comparison of RMSD before and after correction of ft/ 7 4 srt02run01 . 





Root Mean Square Deviation (RMSD) 


RMSD Reduction, % 




Before Correction 


After Correction 




Airspeed (V ), m/s 


2.9855 


0.3373 


88.70 


AoA (a), deg. 


1.3855 


0.3490 


74.81 


AoS(/J), deg. 


2.3143 


0.7520 


67.51 


Roll Angle ( <t> ), deg. 


1.0314 


0.3086 


70.08 


Pitch Angle (6), deg. 


0.5413 


0.2121 


60.82 


Yaw Angle ( y/ ), deg. 


0.8486 


0.4190 


50.62 


Table 2. Comparison of RMSD before and after correction of ft7 4 srt02run0 5. 




Root Mean Square Deviation (RMSD) 


RMSD Reduction, % 




Before Correction 


After Correction 




Airspeed (V ), m/s 


3.3453 


0.0751 


97.76 


AoA ( a ), deg. 


1.7380 


0.4319 


75.15 


AoS(/5), deg. 


2.6914 


0.9741 


63.81 


Roll Angle (^), deg. 


0.7183 


0.3625 


49.53 


Pitch Angle (6>), deg. 


0.4090 


0.2138 


47.73 


Yaw Angle ( y/ ), deg. 


1.6547 


0.3729 


77.46 



9. Conclusions 

This paper presented an easy and straightforward procedure for implementing sensor compatibility 
corrections for fixed- wing UAVs. The procedure implemented the Wiener Type Filter, the RTS 
Smoother, and an airdata noise variances estimation method. The Wiener Type Filter was introduced 
to filter out the high frequency noises of the flight data, thus simplifying the measurement noise matrix 
and process noise matrix estimation. The RTS Smoother serves as the main correction algorithm. It is a 
two pass algorithm that combines the EKF and a backward recursion. The airdata noise variances 
estimation method accomplishes its purpose using the GPS speed measurements. 

The sensor data compatibility correction procedure has been successfully implemented on real 
flight data recorded on the SP-80 UAV of the RMRL. Initial data reconstruction shows that the 
compatibility is poor and the airdata is far much noisier than the AHRS data due to the exposure of the 
airdata probe to the atmospheric turbulence. It was assumed that the AHRS data noises only consist of 
high frequency noises and the noises can be effectively filtered by the Wiener Type Filter. Also, due to 
the high compatibility of the attitude angles, it was assumed that the biases of the angular rates are 
negligible. Since the airdata is corrupted by atmospheric turbulence, the Wiener Type Filter was 
unable to suppress all the noises. Thus, relatively accurate airdata noise variances are needed. The 
airdata noise variances estimation method brings about the concept of data regionality. Under this 
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concept, the data around a particular data point are assumed to have statistical relation with the data 
point. This has allowed the estimation of the expected value and variance at any particular data point 
by taking the data around the data point into account. This concept might not apply universally to all 
problems, but it has adequately fulfilled its purpose for the work presented in this paper. 

After the correction, significant sensor compatibility improvement, especially on the airspeed data 
has been demonstrated. The improvement of the airdata compatibilities is mainly contributed by the 
estimation of the acceleration biases by the RTS. On the other hand, the Wiener Type Filter 
contributes to the improvement of attitude angle compatibilities. 

The main benefit of the sensor data compatibility correction is to provide better data reliability, 
which leads to a better understanding of the dynamics of an air vehicle. Such a correction procedure is 
essential if an accurate dynamic model is to be identified using the flight data. The Wiener Type Filter 
and RTS are both batch (offline) methods. Future improvement to the procedure could include the 
implementation of a real-time (online) algorithm. Being able to perform sensor compatibility 
correction in real-time could help improve the robustness of the autopilot system for both manned and 
unmanned vehicles. 
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